C*********************************************************************C
C* *C
C* chaplit.for *C
C* *C
C* Written by: David L. Huestis, Molecular Physics Laboratory *C
C* *C
C* Copyright (c) 2000 SRI International *C
C* All Rights Reserved *C
C* *C
C* This software is provided on an as is basis; without any *C
C* warranty; without the implied warranty of merchantability or *C
C* fitness for a particular purpose. *C
C* *C
C*********************************************************************C
C*
C* Routines to evaluate formulas from the published literature
C* for the Chapman Function for the effective column depth vs
C* solar zenith angle of an exponential atmosphere.
C*
C* EDIT HISTORY:
C*
C* 01/28/2000 DLH "Best Case" Ch31b evaluation
C*
C* REFERENCES:
C*
C* GLG64 A.E.S. Green, C. S. Lindenmeyer, and M. Griggs,
C* "Molecular Absorption in Planetary Atmospheres,"
C* J. Geophys. Res. _69_, 493-504 (1964).
C*
C* See chapman.for for additional references.
C*
C*********************************************************************C
C ====================================================================
C
C "Best Case" evaluation of Ch(X,chi0) using atm8_chap_Ch31b_22
C for small chi0 and atm8_chap_Ch31b_38 for large chi0
C
C ====================================================================
real*8 function atm8_chap_Ch31b( X, chi0 )
implicit real*8(a-h,o-z)
parameter (rad=57.2957795130823208768d0)
if( (X .le. 0) .or. (chi0 .le. 0) .or. (chi0 .ge. 180) ) then
atm8_chap_Ch31b = 1
return
end if
if( chi0 .gt. 90 ) then
chi = 180 - chi0
else
chi = chi0
end if
sinchi = sin(chi/rad)
coschi = sin((90-chi)/rad)
if( coschi .le. 0 ) then
atm8_chap_Ch31b = atm8_chap_xK1( X )
else if (sinchi .le. 0 ) then
atm8_chap_Ch31b = 1
else if( (2.5d0*log(sinchi) + 4.65d0) .le.
* (log(X) + 5.5d0*log(coschi)) ) then
atm8_chap_Ch31b = atm8_chap_Ch31b_22(X,chi)
else
atm8_chap_Ch31b = atm8_chap_Ch31b_38(X,chi)
end if
if( chi0 .gt. 90 ) then
atm8_chap_Ch31b = 2*exp(X*2*sin((90-chi)/(2*rad))**2)
* * atm8_chap_xK1(X*sinchi) - atm8_chap_Ch31b
end if
return
end
C ====================================================================
C
C Small chi, sec(chi) expansion of Champman [Ch31b, Equation (22)]
C
C Ch(X,chi) = sec(chi) + sum{n=1,5} [ b(n,chi)/X**n ]
C
C which we rewrite as
C
C Ch(X,chi) = sum{n=0.5} [ c(n)/z**n ] * sec(chi)
C
C where
C
C z = X * cos(chi)**2
C c(0) = 1
C c(n) = b(n,chi) * cos(chi)**(2*n), for n > 0
C
C ====================================================================
real*8 function atm8_chap_Ch31b_22( X, chi )
implicit real*8(a-h,o-z)
parameter (rad=57.2957795130823208768d0)
dimension c(0:5)
coschi = cos(chi/rad)
if( (X .le. 0) .or. (chi .le. 0) .or. (chi .ge. 90)
* .or. (coschi .le. 0) ) then
atm8_chap_Ch31b_22 = 1
return
end if
s2 = sin(chi/rad)**2
s4 = s2*s2
c2 = coschi**2
c4 = c2*c2
z = x * c2
c(0) = 1
c(1) = -s2
c(2) = 3*s2
c(3) = -(15*s2+12*c2)*s2
c(4) = (105*s2+60*c2)*s2
c(5) = -(945*s4+1260*s2*c2+360*c4)*s2
sum = c(5)
do i=4,0,-1
sum = sum/z + c(i)
end do
atm8_chap_Ch31b_22 = sum/coschi
return
end
C ====================================================================
C
C Large chi expansion of Champman [Ch31b, Equation (38)]
C
C We write
C
C Ch(X,chi) = y*exp(X)*K1(y)*erfc(sqrt(2*y*V**2))
C + 0.5 * sum{m=0,3} [ Sn(m)/(4*y)**m ]
C
C where
C
C y = X * sin(chi)
C
C V = sqrt( 0.5 * (1/sin(chi) - 1) )
C
C Sn(m) = sum{n=m+1,m+4} [ (2*n-1)...(2*n-2*m+1)
C c(n) * V**(2*n-2*m-1) ]
C
C except that for V .gt. 0.1 we analytically sum the series,
C which significantly improves accuracy for small chi.
C
C Sn(0) = [ (1+2*V**2)/sqrt(1+V**2) -1 ] / V
C
C Sn(1) = [ (d/dV) Sn(0) - c(1) ] / V
C
C Sn(2) = [ (d/dV) Sn(1) - 3*c(2) ] / V
C
C Sn(3) = [ (d/dV) Sn(2) - 15*c(3) ] / V
C
C With this formulations, the Sn(m) vanish as V for small
C V, and as 1/V for large V.
C
C ====================================================================
real*8 function atm8_chap_Ch31b_38( X, chi )
implicit real*8(a-h,o-z)
parameter (rad=57.2957795130823208768d0)
parameter (small=1.0d-2)
dimension c(7), Sn(0:3)
common/atm8_chap_cm/Fn(0:3)
data c/1.5d0, -6.25d-1, 4.375d-1, -3.515625d-1,
* 3.0078125d-1, -2.666015625d-1, 2.4169921875d-1/
if( (X .le. 0) .or. (chi .le. 0) .or. (chi .gt. 90) ) then
atm8_chap_Ch31b_38 = 1
return
end if
sinchi = sin(chi/rad)
y = X * sinchi
beta = X * 2*sin((90-chi)/(2*rad))**2
C beta = X*(1-sinchi) = 2 * y * V**2
Fn(0) = atm8_chap_xK1( y ) * atm8_chap_xerfc( sqrt(beta) )
V2 = 0.5d0 * beta / y
if( V2 .le. 0 ) then
do n=1,3
Fn(n) = Fn(n-1)
end do
else
V = sqrt( V2 )
if( V2 .le. small ) then
Sn(0) = V*(c(1) + V2*(c(2) + V2*(c(3) + V2*c(4))))
Sn(1) = V*(3*c(2) + V2*(5*c(3) + V2*(7*c(4) + V2*9*c(5))))
Sn(2) = V*(15*c(3) + V2*(35*c(4) + V2*(63*c(5)
* + V2*99*c(6))))
Sn(3) = V*(105*c(4) + V2*(315*c(5) + V2*(693*c(6)
* + V2*1287*c(7))))
else
V12 = 1+V2
V3 = V2*V
V4 = V3*V
V5 = V4*V
V6 = V5*V
r = sqrt(V12)
r3 = r*V12
r5 = r3*V12
r7 = r5*V12
C h(V) = (1+2*V**2)/sqrt(1+V**2) - 1
h = ( 2*r - 1/r -1 )
hp = V*(3+2*V2)/r3
hpp = 3/r5
hppp = -15*V/r7
Sn(0) = h / V
Sn(1) = ( -h/V2 + hp/V - c(1) )/V
Sn(2) = ( +3*h/V4 -3*hp/V3 +hpp/V2 +c(1)/V2 -3*c(2) )/V
Sn(3) = ( -15*h/V6 +15*hp/V5 -6*hpp/V4 + hppp/V3
* -3*c(1)/V4 +3*c(2)/V2 -15*c(3) )/V
end if
Y4 = 4*Y
Y42 = Y4*Y4
Y43 = Y42*Y4
Fn(0) = Fn(0) + 0.5d0 * Sn(0)
Fn(1) = Fn(0) + 0.5d0 * Sn(1) / Y4
Fn(2) = Fn(1) + 0.5d0 * Sn(2) / Y42
Fn(3) = Fn(2) + 0.5d0 * Sn(3) / Y43
end if
atm8_chap_Ch31b_38 = Fn(3)
return
end
C ====================================================================
C
C Empirical Equation (3) of Fitmaurice [Fi64]
C
C ====================================================================
real*8 function atm8_chap_Fi64_3( X, theta )
implicit real*8(a-h,o-z)
parameter (pi = 3.1415926535897932384d0)
parameter (rad=57.2957795130823208768d0)
atm8_chap_Fi64_3 = sqrt(X*pi/2)
* * atm8_chap_xerfc( sqrt(X/2)*cos(theta/rad) )
return
end
C ====================================================================
C
C Empirical sec(chi) expansion, Equation (5) of Fitmaurice [Fi64]
C
C ====================================================================
real*8 function atm8_chap_Fi64_5( X, theta )
implicit real*8(a-h,o-z)
parameter (rad=57.2957795130823208768d0)
cost = cos(theta/rad)
if( (X .le. 0) .or. (theta .le. 0) .or. (theta .gt. 90)
* .or. (cost .le. 0) ) then
atm8_chap_F64_5 = 1
return
end if
t = 1/( X * cost**2 )
atm8_chap_Fi64_5 = (1 + t*( -1 + t*( 3 - t*15) ) ) / cost
return
end
C ====================================================================
C
C Empirical Equation (11) of Green et al [GLG64]
C
C ====================================================================
real*8 function atm8_chap_GLG64_11( X, theta )
implicit real*8(a-h,o-z)
parameter (pi = 3.1415926535897932384d0)
parameter (pi2 = pi/2)
parameter (rad=57.2957795130823208768d0)
if( (X .le. 0) .or. (theta .le. 0) .or. (theta .gt. 90) ) then
atm8_chap_GLG64_11 = 1
return
end if
alpha = ( 1 - 0.115d0*pi2**2
* - 0.5d0*(pi2**2)/dlog(sqrt(pi2*X)) ) / (pi2**4)
theta2 = (theta/rad)**2
bot = 1+theta2*(-0.115d0-theta2*alpha)
if( bot .gt. 0 ) then
atm8_chap_GLG64_11 = exp( 0.5d0*theta2/bot )
else
atm8_chap_GLG64_11 = 1
end if
return
end
C ====================================================================
C
C Empirical Equation (52) of Swider [Sw64], Corrected [SG69]
C
C ====================================================================
real*8 function atm8_chap_Sw64_52( X, chi0 )
implicit real*8(a-h,o-z)
parameter (pi = 3.1415926535897932384d0)
parameter (rad=57.2957795130823208768d0)
sinchi = sin(chi0/rad)
t0sqr = 2*sin( (90-chi0)/(2*rad) )**2
XC = sqrt(X*t0sqr)
atm8_chap_Sw64_52 = -X*cos(chi0/rad) + sqrt(1+X*(1+sinchi))
* *( XC + (sqrt(pi)/2)*atm8_chap_xerfc(XC) )
return
end